library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(tidyr)
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
devtools::load_all('..')
## ℹ Loading hector
data_dir <- 'data/'
ini_dir <- 'data/pic_hector_ini/climate/'
This is true NBP: atmosphere to land, positive NBP = a sink
# load the file passed from Dawn
gcp_data_dawn <- read.csv(paste0(data_dir, "nbp_gcp.csv"))
# Load the file pulled directly from GCP:
read.csv(paste0(data_dir,'Global_Carbon_Budget_2022v1p0_historical_co2_tab.csv'), stringsAsFactors = F) %>%
select(year=Year, e_luc=land.use.change.emissions, s_land = land.sink, Units) %>%
tidyr::replace_na(list(e_luc=0)) %>%
# If we want to sum e_luc + s_land, uncomment from the `gather` to `ungroup`
# calls here and comment out the `mutate`:
# gather(reporting_id, value, -year, -Units) %>%
# group_by(year, Units) %>%
# summarize(nbp = -sum(value, na.rm=T)) %>%
# ungroup %>%
# Following slide 56 of GCP_data/papers/GCP_CarbonBudget_2023_slides_v1.0-2-Alissa-Haward.pdf
mutate(nbp = (s_land - e_luc)) ->
gcp_data
Check if nbp = (s_land-e_luc) from the raw GCP download
agrees with file from Dawn. Note Dawn’s file did
-(s_land-e_luc) to have things in the land to atmosphere
direction instead, and compare with offline model outputs directly. so
have to do minus a minus
print("Difference to Dawn's file:")
## [1] "Difference to Dawn's file:"
print(max(abs(gcp_data$nbp - -gcp_data_dawn$nbp)))
## [1] 0.01
Label data and convert units
# Label data and convert units
gcp_data$scenario <- "1. Global Carbon Project"
gcp_data$nbp_raw <- gcp_data$nbp*1000
gcp_data$nbp <- zoo::rollmean(gcp_data$nbp_raw,k=10, align='right', fill= NA)
gcp_data %>%
filter(year >=1970)
## year e_luc s_land Units nbp scenario nbp_raw
## 1 1970 1.26 0.71 GtC/yr -228 1. Global Carbon Project -550
## 2 1971 1.25 2.45 GtC/yr -25 1. Global Carbon Project 1200
## 3 1972 1.25 1.25 GtC/yr 4 1. Global Carbon Project 0
## 4 1973 1.21 1.89 GtC/yr 139 1. Global Carbon Project 680
## 5 1974 1.18 4.19 GtC/yr 420 1. Global Carbon Project 3010
## 6 1975 1.19 2.65 GtC/yr 670 1. Global Carbon Project 1460
## 7 1976 1.17 2.99 GtC/yr 844 1. Global Carbon Project 1820
## 8 1977 1.20 1.51 GtC/yr 841 1. Global Carbon Project 310
## 9 1978 1.15 2.78 GtC/yr 900 1. Global Carbon Project 1630
## 10 1979 1.11 1.42 GtC/yr 987 1. Global Carbon Project 310
## 11 1980 1.11 0.41 GtC/yr 972 1. Global Carbon Project -700
## 12 1981 1.20 2.22 GtC/yr 954 1. Global Carbon Project 1020
## 13 1982 1.20 1.67 GtC/yr 1001 1. Global Carbon Project 470
## 14 1983 1.31 0.28 GtC/yr 830 1. Global Carbon Project -1030
## 15 1984 1.45 2.89 GtC/yr 673 1. Global Carbon Project 1440
## 16 1985 1.39 2.77 GtC/yr 665 1. Global Carbon Project 1380
## 17 1986 1.40 2.31 GtC/yr 574 1. Global Carbon Project 910
## 18 1987 1.40 0.35 GtC/yr 438 1. Global Carbon Project -1050
## 19 1988 1.38 2.06 GtC/yr 343 1. Global Carbon Project 680
## 20 1989 1.33 3.57 GtC/yr 536 1. Global Carbon Project 2240
## 21 1990 1.33 2.33 GtC/yr 706 1. Global Carbon Project 1000
## 22 1991 1.30 2.25 GtC/yr 699 1. Global Carbon Project 950
## 23 1992 1.33 2.26 GtC/yr 745 1. Global Carbon Project 930
## 24 1993 1.33 2.86 GtC/yr 1001 1. Global Carbon Project 1530
## 25 1994 1.45 1.47 GtC/yr 859 1. Global Carbon Project 20
## 26 1995 1.44 1.77 GtC/yr 754 1. Global Carbon Project 330
## 27 1996 1.44 3.40 GtC/yr 859 1. Global Carbon Project 1960
## 28 1997 1.94 3.19 GtC/yr 1089 1. Global Carbon Project 1250
## 29 1998 1.50 1.51 GtC/yr 1022 1. Global Carbon Project 10
## 30 1999 1.50 3.59 GtC/yr 1007 1. Global Carbon Project 2090
## 31 2000 1.40 3.75 GtC/yr 1142 1. Global Carbon Project 2350
## 32 2001 1.29 2.38 GtC/yr 1156 1. Global Carbon Project 1090
## 33 2002 1.41 0.99 GtC/yr 1021 1. Global Carbon Project -420
## 34 2003 1.54 2.31 GtC/yr 945 1. Global Carbon Project 770
## 35 2004 1.46 3.47 GtC/yr 1144 1. Global Carbon Project 2010
## 36 2005 1.29 1.88 GtC/yr 1170 1. Global Carbon Project 590
## 37 2006 1.38 3.11 GtC/yr 1147 1. Global Carbon Project 1730
## 38 2007 1.21 2.71 GtC/yr 1172 1. Global Carbon Project 1500
## 39 2008 1.27 3.63 GtC/yr 1407 1. Global Carbon Project 2360
## 40 2009 1.37 2.99 GtC/yr 1360 1. Global Carbon Project 1620
## 41 2010 1.32 3.32 GtC/yr 1325 1. Global Carbon Project 2000
## 42 2011 1.35 4.11 GtC/yr 1492 1. Global Carbon Project 2760
## 43 2012 1.32 2.32 GtC/yr 1634 1. Global Carbon Project 1000
## 44 2013 1.26 3.56 GtC/yr 1787 1. Global Carbon Project 2300
## 45 2014 1.34 3.66 GtC/yr 1818 1. Global Carbon Project 2320
## 46 2015 1.47 2.01 GtC/yr 1813 1. Global Carbon Project 540
## 47 2016 1.24 2.69 GtC/yr 1785 1. Global Carbon Project 1450
## 48 2017 1.18 3.56 GtC/yr 1873 1. Global Carbon Project 2380
## 49 2018 1.14 3.65 GtC/yr 1888 1. Global Carbon Project 2510
## 50 2019 1.24 3.04 GtC/yr 1906 1. Global Carbon Project 1800
## 51 2020 1.11 3.11 GtC/yr 1906 1. Global Carbon Project 2000
## 52 2021 1.08 3.45 GtC/yr 1867 1. Global Carbon Project 2370
#comparison with GCP
# package scenario
rcp <- "ssp245"
scenario_file <- paste0("input/hector_",rcp,".ini")
ini_file245 <- system.file(scenario_file, package="hector")
#### QUESTION - is this right hector package if did devtools::load_all
#### ANSWER: yes, because added different ini to inst/input directory and after
#### devtools::load_all(..) that new ini is then accessible.
core245 <- hector::newcore(ini_file245, suppresslogging = FALSE, name='245')
hector::run(core245,runtodate = 2100)
## Hector core: 245
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini
out <- hector::fetchvars(core245, dates=1746:2100, vars=c(hector::NBP()))
gcam_emiss <- read.csv('/Users/snyd535/Documents/GitHub/hector/inst/input/gcam_emissions.csv', skip =4)
ggplot(gcam_emiss) +
geom_line(mapping = aes(x=Date, y=luc_emissions))
library(rgcam)
get_gcam_luc_emissR <- function(db_name="database_basexdb", gcam_dir, scenario,
read_from_file=FALSE,
filename="data/gcam_land_alloc.csv"){
if (read_from_file) {
gcam_land_alloc <- read.csv2(file=filename,header=TRUE)
}
else {
base_conn <- rgcam::localDBConn(gcam_dir, db_name)
land_alloc_query <- '<query title="LUC emissions by region">
<axis1 name="LandLeaf">LandLeaf</axis1>
<axis2 name="Year">land-use-change-emission[@year]</axis2>
<xPath buildList="true" dataName="land-use-change-emission" group="false" sumAll="true">/LandNode[@name=\'root\' or @type=\'LandNode\' (:collapse:)]//
land-use-change-emission[@year>1970]/text()</xPath>
<comments/>
</query>'
new.proj <- rgcam::addSingleQuery(base_conn, "new.proj", "LUC_R", land_alloc_query, c(scenario), clobber=TRUE)
gcam_land_alloc <- rgcam::getQuery(new.proj, "LUC_R")
# write.csv(gcam_land_alloc, file=filename, row.names = FALSE)
}
return(gcam_land_alloc)
}
luc_R<- get_gcam_luc_emissR(db_name = 'database_basexdbGCAM',
gcam_dir = 'data/pic_hector_ini',
scenario = 'GCAM')
## Database scenarios: GCAM
ggplot(luc_R %>% filter(landleaf == 'Pasture_MissouriR', region == 'USA', year <=2015)) +
geom_line(mapping = aes(x=year, y=value, color = scenario)) +
ylab('MtC/yr')+
ggtitle('Pasture_MissouriR GCAM output LUC emissions')
ggplot(luc_R %>%
group_by(Units, scenario, year) %>%
summarize(value = sum( value)) %>%
ungroup %>%
filter(year<=2015)) +
geom_point(mapping = aes(x=year, y=value, color = scenario)) +
ylab('MtC/yr')+
ggtitle('Global aggregate GCAM output LUC emissions') +
geom_line(data = gcam_emiss %>% filter(Date >= 1975, Date <= 2015),
mapping = aes(x=Date, y=luc_emissions*1e3), color = 'black')
## `summarise()` has grouped output by 'Units', 'scenario'. You can override using
## the `.groups` argument.
gcp_data %>%
select(year, scenario, nbp, Units) %>%
select(year, value=nbp, units=Units) %>%
mutate(variable = 'NBP',
scenario ="GCP",
value = 0.001*value,
units = 'PgC/year'
) %>%
bind_rows(out) ->
out2
ggplot(out2 %>% filter(year <= 2021)) +
geom_line(mapping = aes(x=year,y=value, colour=scenario)) +
geom_line(data= (gcp_data %>%
select(year, scenario, nbp_raw, Units) %>%
select(year, value=nbp_raw, units=Units) %>%
mutate(variable = 'NBP_raw',
scenario ="GCP_raw",
value = 0.001*value,
units = 'PgC/year'
)),
mapping = aes(x=year, y=value), color = 'grey', alpha = 0.5) +
ggtitle('gray is raw GCP')
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_line()`).
rm(out)
rm(out2)
rm(gcp_data_dawn)
rm(luc_R)
The 0-level test for whether the NBP constraint works is if you can do a Hector run, extract the NBP from that run, use that NBP as a constraint for a new Hector run, and get the same results.
This is trickier than it sounds, due to uninformative errors
triggered in fluxpool.hpp:
inline void fluxpool::set(double v, unit_types u, bool track = false,
string pool_name = "test") {
name = pool_name;
if (v < 0) {
H_ASSERT(v >= 0, "Flux and pool values may not be negative in 2 " + name);
}
Solution Part 1: only use an NBP constraint with positive values.
In practice, this just means only using years >1970 or so as the constraint.
for GCP, depending on how you smooth it, this can be 1972. But if you try to use 1972 on with hector data, need to figure out a different diffusivity. very dumb.
This led to some different errors, sometimes, that made me think I was on the right track but still ultimately gave the same flux values may not be negative error most of the time.
Solution Part 2: Lower the ocean heat diffusivity from its default value of 2.38 cm2/s to 2.31 cm2/s (for a 245).
ACS had run some NBP constraints a year+ ago and they were NOT as finicky to get to solve, but that version of hector had different beta, q10, and diffusivity.
the old (1.16) vs new value of diffusivity was the biggest difference (new was about double old), so started at 2.38 and just started incremementing lower until the constraint solved.
definitely had some runs with 2.34 that solved but can’t
replicate. I think this is not a real instability in hector
code but is instead user error - in rstudio, you can’t cmd+enter a bunch
of rhector function calls together - some end up kind of
skipped because the newcore and setvar
functions seem to take an extra moment and need to be run on their own
before proceeding?
but also only sometimes. Seems to always be fine when knitting at least.
handful of steps with the above, so make a function instead of copying and pasting every time.
adjusting ocean heat diffusivity to 2.31 cm2/s
adjust_hector_run <- function(initial_ini,
scen_name,
diffusivity = 2.31){
rcp <- initial_ini
scenario_file <- paste0("input/hector_",rcp,".ini")
ini_file <- system.file(scenario_file, package="hector")
core <- hector::newcore(ini_file, suppresslogging = FALSE, name=scen_name)
hector::run(core,runtodate = 2100)
setvar(core, NA, DIFFUSIVITY(), diffusivity, "cm2/s") # default 2.38
print(fetchvars(core, NA, DIFFUSIVITY()))
hector::run(core,runtodate = 2100)
return(core)
}
output_extractor <- function(core){
fetchvars(core,
dates=1750:2100,
vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE(),
SOIL_C(),
VEG_C())) %>%
mutate(variable = if_else( variable == 'NBP', 'A.NBP',
if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
if_else(variable == 'global_tas', 'C.global_tas',
if_else(variable == 'LL_ocean_uptake', 'F.LL_ocean_uptake',
if_else(variable == 'HL_ocean_uptake', 'G.HL_ocean_uptake',
if_else( variable == 'soil_c', 'D.soil_c',
if_else(variable == 'veg_c', 'E.veg_c',
variable)))))))) %>%
mutate(variable = paste0(variable, '~', units))
}
# Default hector
rcp <- "ssp245"
scenario_file <- paste0("input/hector_",rcp,".ini")
ini_file245 <- system.file(scenario_file, package="hector")
core245 <- hector::newcore(ini_file245, name = 'ssp245_default')
hector::run(core245,runtodate = 2100)
## Hector core: ssp245_default
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini
# Core with adjusted params that work for constraint
#reset(core_test)
core_test <- adjust_hector_run(initial_ini = 'ssp245',
scen_name = 'ssp245_diff2.31')
## scenario year variable value units
## 1 ssp245_diff2.31 NA diff 2.31 cm2/s
## Auto-resetting core to 0
hector::run(core_test,runtodate = 2100)
## Hector core: ssp245_diff2.31
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini
# extract outputs
bind_rows(output_extractor(core245),
output_extractor(core_test)) ->
plot_out
ggplot(plot_out) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
ggtitle('Effect of diffusivity 2.31 cm2/s (old = 2.38)')
plot_out %>%
filter(scenario == 'ssp245_default') %>%
rename(ssp245_default = value) %>%
select(-scenario) %>%
left_join(plot_out %>%
filter(scenario == 'ssp245_diff2.31') %>%
rename(ssp245_diff2.31 = value) %>%
select(-scenario),
by = c('year', 'variable', 'units') ) %>%
group_by(variable, units) %>%
summarize(maxdiff = max(abs(ssp245_default-ssp245_diff2.31)),
maxnormdiff = max(abs(ssp245_default-ssp245_diff2.31)) /max(abs(ssp245_default))) %>%
ungroup %>%
knitr::kable()
## `summarise()` has grouped output by 'variable'. You can override using the
## `.groups` argument.
| variable | units | maxdiff | maxnormdiff |
|---|---|---|---|
| A.NBP~Pg C/yr | Pg C/yr | 0.0076322 | 0.0024744 |
| B.CO2_concentration~ppmv CO2 | ppmv CO2 | 0.3029845 | 0.0005408 |
| C.global_tas~degC | degC | 0.0132185 | 0.0048332 |
| D.soil_c~Pg C | Pg C | 0.4119093 | 0.0001965 |
| E.veg_c~Pg C | Pg C | 0.1063067 | 0.0001582 |
| F.LL_ocean_uptake~Pg C | Pg C | 0.0025119 | 0.0009716 |
| G.HL_ocean_uptake~Pg C | Pg C | 0.0005476 | 0.0002809 |
ggplot(plot_out %>% filter(grepl('ocean', variable))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
ggtitle('Effect of diffusivity 2.31 cm2/s (old = 2.38)')
So not surprisingly based on a lot of other things we have seen over the last couple years with hector, the ocean is the part that needs maintenance.
good to know diffusivity is a lever. Will do some sensitivity tests further down sweeping out different diffusivity values vs different magnitudes of NBP (will explain more there).
zoom in some more on terrestrial carbon quantities
ggplot(plot_out %>% filter(grepl('NBP', variable) | grepl('veg', variable) | grepl('soil', variable)) )+
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=3) +
ggtitle('Effect of diffusivity 2.31 cm2/s (old = 2.38)')
So from this point, all Hector runs will be done with 2.31cm2/s for diffusivity so the NBP constraints will actually run without running into negative flux pools.
If we extract the entire nbp trajectory from a hector run and make it the constraint for a run with otherwise identical settings, do we get identical results?
# clean up some variables
rm(core_test)
rm(core245)
rm(plot_out)
# initialize a holder of outputs
output_holder <- list()
# Do the run to generate the constraing
core <- adjust_hector_run(initial_ini = 'ssp245',
scen_name = 'ssp245_gen')
## scenario year variable value units
## 1 ssp245_gen NA diff 2.31 cm2/s
## Auto-resetting core to 0
output_holder[['ssp245_gen']] <- output_extractor(core)
rm(core)
# Extract Constraint
output_holder[['ssp245_gen']] %>%
filter(grepl('NBP', variable),
year >= 1970,
year<= 2100) ->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp245',
scen_name = 'ssp245_1970~2100')
## scenario year variable value units
## 1 ssp245_1970~2100 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1969
## Hector core: ssp245_1970~2100
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini
output_holder[['ssp245_1970~2100']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
ggtitle('testing NBP constraint SSP245')
Previously only compared NBP before and after applying constraint.
clear from looking at these additional quantities, it’s not working correctly.
differences in soil and veg carbon are possibly a bit weird (depending on how constraint implemented in code). But in theory, calculations from the land box don’t go anywhere when there is an NBP constraint, so a difference is ok as long as the land calculations don’t co anywhere.
But the ocean uptakes differing as well suggests the land calculations ARE going somewhere and being used instead of the constraint.
HYPOTHESIS in code, the NBP constraint is overwriting some values in some places (which is why NBP constraint returns identical) but probably not in all of the places it should.
So constraint being on = land box evolves differently, with lower soil and veg carbon stocks. Lower land carbon stocks (soil especially) means lower ocean uptake in reaction? is that right directions?
## `summarise()` has grouped output by 'variable'. You can override using the
## `.groups` argument.
| variable | units | maxdiff | maxnormdiff |
|---|---|---|---|
| A.NBP~Pg C/yr | Pg C/yr | 0.0000000 | 0.0000000 |
| B.CO2_concentration~ppmv CO2 | ppmv CO2 | 4.2906869 | 0.0076546 |
| C.global_tas~degC | degC | 0.0209422 | 0.0076212 |
| D.soil_c~Pg C | Pg C | 52.2774744 | 0.0249493 |
| E.veg_c~Pg C | Pg C | 8.2995916 | 0.0123482 |
| F.LL_ocean_uptake~Pg C | Pg C | 0.1860536 | 0.0719548 |
| G.HL_ocean_uptake~Pg C | Pg C | 0.1499027 | 0.0768922 |
Know ssp245 can be weird (github issue), so double check other scenarios
# Do the run to generate the constraing
core <- adjust_hector_run(initial_ini = 'ssp126',
scen_name = 'ssp126_gen')
## scenario year variable value units
## 1 ssp126_gen NA diff 2.31 cm2/s
## Auto-resetting core to 0
output_holder[['ssp126_gen']] <- output_extractor(core)
rm(core)
# Extract Constraint
output_holder[['ssp126_gen']] %>%
filter(grepl('NBP', variable),
year >= 1970,
year<= 2100) ->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp126',
scen_name = 'ssp126_1970~2100')
## scenario year variable value units
## 1 ssp126_1970~2100 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1969
## Hector core: ssp126_1970~2100
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp126.ini
output_holder[['ssp126_1970~2100']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('126', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
ggtitle('testing NBP constraint SSP126')
# Do the run to generate the constraing
core <- adjust_hector_run(initial_ini = 'ssp585',
scen_name = 'ssp585_gen')
## scenario year variable value units
## 1 ssp585_gen NA diff 2.31 cm2/s
## Auto-resetting core to 0
output_holder[['ssp585_gen']] <- output_extractor(core)
rm(core)
# Extract Constraint
output_holder[['ssp585_gen']] %>%
filter(grepl('NBP', variable),
year >= 1970,
year<= 2100) ->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
scen_name = 'ssp585_1970~2100')
## scenario year variable value units
## 1 ssp585_1970~2100 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1969
## Hector core: ssp585_1970~2100
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini
output_holder[['ssp585_1970~2100']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
ggtitle('testing NBP constraint SSP585')
Seems useful to understand what happens when you just constrain the (usable) historical period.
usable in the sense that need NBP constraint values not to be negative, so years > 1970.
cutting at 2021 because that’s the last year in the GCP data I have.
# Constraint with updated params
output_holder[['ssp245_gen']] %>%
filter(grepl('NBP', variable),
year >= 1970,
year<= 2021) ->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp245',
scen_name = 'ssp245_1970~2021')
## scenario year variable value units
## 1 ssp245_1970~2021 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1969
## Hector core: ssp245_1970~2021
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini
output_holder[['ssp245_1970~2021']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
ggtitle('testing NBP constraint SSP245')
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
By terminating the constraint in 2021, all variables are able to drift closer to their default hector behavior.
NOTE the abrupt jump in NBP when the constraint ceases. HYPOTHESIS The background land calculations going on while the constraint is in place are immediately what kicks in at the end of the constraint, rather than the constraint value.
In 2020, background calculation produces lower carbon stocks in land = abrupt jump in land sink? can’t tell if that’s right but blue vs green and red vs green curve comparison in this figure is the key to figuring it out, I think.
A few years in the GCP data from 1970-1980 have nbp<0.
This seems to cause the NBP constraint Hector actually sees to remain at 0, leading to a huge rebound spike when the constraint ends in 2021.
Just droping the specific years with negative nbp causes Hector to fill in 0s for those years, leading to a rebound at each of these years (not shown).
# Run with GCP NBP
gcp_data %>%
filter(year >= 1972,
year<= 2021) %>%
#filter(nbp>0) %>%
mutate(value = 0.001*nbp,
Units = "Pg C/yr")->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp245',
scen_name = 'ssp245_GCP1972~2021')
## scenario year variable value units
## 1 ssp245_GCP1972~2021 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1971
## Hector core: ssp245_GCP1972~2021
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp245.ini
output_holder[['ssp245_GCP1972~2021']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
ggtitle('testing NBP constraint SSP245')
ggplot(do.call(rbind, output_holder) %>% filter(grepl('245', scenario), year >= 1965, year <=2050)) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
ggtitle('testing NBP constraint SSP245')
# Constraint with updated params
output_holder[['ssp585_gen']] %>%
filter(grepl('NBP', variable),
year >= 1970,
year<= 2021) ->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
scen_name = 'ssp585_1970~2021')
## scenario year variable value units
## 1 ssp585_1970~2021 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1969
## Hector core: ssp585_1970~2021
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini
output_holder[['ssp585_1970~2021']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
ggtitle('testing NBP constraint SSP585')
# Run with GCP NBP
gcp_data %>%
filter(year >= 1972,
year<= 2021) %>%
#filter(nbp>0) %>%
mutate(value = 0.001*nbp,
Units = "Pg C/yr")->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
scen_name = 'ssp585_GCP1972~2021')
## scenario year variable value units
## 1 ssp585_GCP1972~2021 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1971
## Hector core: ssp585_GCP1972~2021
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini
output_holder[['ssp585_GCP1972~2021']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
ggtitle('testing NBP constraint SSP585')
# Run with GCP NBP
gcp_data %>%
filter(year >= 1972,
year<= 2000) %>%
#filter(nbp>0) %>%
mutate(value = 0.001*nbp,
Units = "Pg C/yr")->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp585',
scen_name = 'ssp585_GCP1972~2000')
## scenario year variable value units
## 1 ssp585_GCP1972~2000 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1971
## Hector core: ssp585_GCP1972~2000
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp585.ini
output_holder[['ssp585_GCP1972~2000']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario))) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
geom_vline(xintercept = 2021, size = 0.1, color = 'grey')+
ggtitle('testing NBP constraint SSP585')
ggplot(do.call(rbind, output_holder) %>% filter(grepl('585', scenario), year >= 1965, year <=2050)) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
geom_vline(xintercept = 2021, size = 0.1, color = 'grey')+
ggtitle('testing NBP constraint SSP585')
# Constraint with updated params
output_holder[['ssp126_gen']] %>%
filter(grepl('NBP', variable),
year >= 1970,
year<= 2021) ->
nbp_constrainer
# Set up core for constraint, set constraint, run wtih constraint
core <- adjust_hector_run(initial_ini = 'ssp126',
scen_name = 'ssp126_1970~2021')
## scenario year variable value units
## 1 ssp126_1970~2021 NA diff 2.31 cm2/s
## Auto-resetting core to 0
setvar(core, nbp_constrainer$year, NBP_CONSTRAIN(), nbp_constrainer$value, "Pg C/yr")
hector::run(core,runtodate = 2100)
## Auto-resetting core to 1969
## Hector core: ssp126_1970~2021
## Start date: 1745
## End date: 2300
## Current date: 2100
## Input file: /Users/snyd535/Documents/GitHub/hector/inst/input/hector_ssp126.ini
output_holder[['ssp126_1970~2021']] <- output_extractor(core)
rm(core)
ggplot(do.call(rbind, output_holder) %>% filter(!grepl('GCP', scenario), year >= 1965, year <=2050)) +
geom_line(mapping = aes(x=year, y = value, color = scenario, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=2) +
geom_vline(xintercept = 2021, size = 0.2)+
ggtitle('testing NBP constraint')
So my guess of what’s happening, just based on the plots, is the land calculations still run in the background and the NBP constraint is supposed to just ignore them. But somewhere, the land calculations aren’t actually being ignored and the ocean uptake reacts to them.
So I think Atm C is getting the nbp constraint for land, but also the ocean differences reacting to background land caluclations insteadof NBP and that’s what makes AtmC a little different?
No idea:
A. Where in the code the land calculations that should be ignored are getting passed on to anything
B. Why those background land calculations are different
There is also:
Can update the below when NBP constraint actually mechanically works.
all of above was just trying to get NBP cosntraint to run at all.
Need:
to make sure that for an ordered set of NBP constraints, T and CO2 should be following sensible order
sensitivity on magnitude of constraint vs diffusivity
keep it simple - couple multiples of 245 NBP with 2.31
# Get constraint:
core_constrain <- hector::newcore(ini_file245, name = '245_diff2.31_generate')
setvar(core_constrain, NA, DIFFUSIVITY(), 2.31, "cm2/s") # default 2.38
hector::run(core_constrain,runtodate = 2050)
nbp <- fetchvars(core_constrain, 1965:2050, NBP())
out <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=1)
# Constraint with updated params
nbp %>%
filter(scenario == '245_diff2.31_generate',
year >= 1970,
year<= 2021) ->
nbp_constrainer
for (k in c(0.5, 1, 2, 5, 8)){
core_k <- hector::newcore(ini_file245, name = paste0('245_diff2.31_constrain_k',k))
setvar(core_k, NA, DIFFUSIVITY(), 2.31, "cm2/s")
setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
hector::run(core_k,runtodate = 2050)
fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=k) %>%
bind_rows(out) ->
out
rm(core_k)
}
# Order the variables for faceting
out %>%
mutate(variable = if_else( variable == 'NBP', 'A.NBP',
if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
if_else(variable == 'global_tas', 'C.global_tas',
if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
variable)))))) %>%
mutate(variable = paste0(variable, '~', units)) ->
out2
ggplot(out2) +
geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=1) +
geom_vline(xintercept = 2021, size = 0.3) +
scale_color_continuous() +
ggtitle('diffusivity 2.31 cm2/s')
keep it simple - couple multiples of 245 NBP with 2.31
# Get constraint:
nbp <- fetchvars(core_constrain, 1965:2050, NBP())
out <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=1)
# Constraint with updated params
nbp %>%
filter(scenario == '245_diff2.31_generate',
year >= 1970,
year<= 2050) ->
nbp_constrainer
for (k in c(0.5, 1, 2, 5)){
core_k <- hector::newcore(ini_file245, name = paste0('245_diff2.31_constrain_k',k))
setvar(core_k, NA, DIFFUSIVITY(), 2.31, "cm2/s")
setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
hector::run(core_k,runtodate = 2050)
fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=k) %>%
bind_rows(out) ->
out
rm(core_k)
}
# Order the variables for faceting
out %>%
mutate(variable = if_else( variable == 'NBP', 'A.NBP',
if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
if_else(variable == 'global_tas', 'C.global_tas',
if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
variable)))))) %>%
mutate(variable = paste0(variable, '~', units)) ->
out2
ggplot(out2) +
geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=1) +
# geom_vline(xintercept = 2021, size = 0.3) +
scale_color_continuous() +
ggtitle('diffusivity 2.31 cm2/s')
use 1.16, show that nbp/co2/temp don’t behave the way they should.
# Get constraint:
core_constrain <- hector::newcore(ini_file245, name = '245_diff1.16_generate')
setvar(core_constrain, NA, DIFFUSIVITY(), 1.16, "cm2/s") # default 2.38
hector::run(core_constrain,runtodate = 2050)
nbp <- fetchvars(core_constrain, 1965:2050, NBP())
out <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=1)
# Constraint with updated params
nbp %>%
filter(scenario == '245_diff1.16_generate',
year >= 1970,
year<= 2021) ->
nbp_constrainer
for (k in c(0.5, 1, 2, 5, 8)){
core_k <- hector::newcore(ini_file245, name = paste0('245_diff1.16_constrain_k',k))
setvar(core_k, NA, DIFFUSIVITY(), 1.16, "cm2/s")
setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
hector::run(core_k,runtodate = 2050)
fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=k) %>%
bind_rows(out) ->
out
rm(core_k)
}
# Order the variables for faceting
out %>%
mutate(variable = if_else( variable == 'NBP', 'A.NBP',
if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
if_else(variable == 'global_tas', 'C.global_tas',
if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
variable)))))) %>%
mutate(variable = paste0(variable, '~', units)) ->
out2
ggplot(out2) +
geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=1) +
geom_vline(xintercept = 2021, size = 0.3) +
scale_color_continuous() +
ggtitle('diffusivity 1.16 cm2/s')
use 1.16, show that nbp/co2/temp don’t behave the way they should.
# Get constraint:
nbp <- fetchvars(core_constrain, 1965:2050, NBP())
out <-fetchvars(core_constrain,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=1)
# Constraint with updated params
nbp %>%
filter(scenario == '245_diff1.16_generate',
year >= 1970,
year<= 2050) ->
nbp_constrainer
for (k in c(0.5, 1, 2, 5)){
core_k <- hector::newcore(ini_file245, name = paste0('245_diff1.16_constrain_k',k))
setvar(core_k, NA, DIFFUSIVITY(), 1.16, "cm2/s")
setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$value, "Pg C/yr")
hector::run(core_k,runtodate = 2050)
fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=k) %>%
bind_rows(out) ->
out
rm(core_k)
}
# Order the variables for faceting
out %>%
mutate(variable = if_else( variable == 'NBP', 'A.NBP',
if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
if_else(variable == 'global_tas', 'C.global_tas',
if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
variable)))))) %>%
mutate(variable = paste0(variable, '~', units)) ->
out2
ggplot(out2) +
geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=1) +
geom_vline(xintercept = 2021, size = 0.3) +
scale_color_continuous() +
ggtitle('diffusivity 1.16 cm2/s')
# GCP
gcp_data %>%
filter(year >= 1972,
year<= 2021) %>%
#filter(nbp>0) %>%
mutate(nbp = 0.001*nbp,
Units = "Pg C/yr")->
nbp_constrainer
coreGCP <- hector::newcore(ini_file245, name = 'GCP_diff2.31_constraint_1972~2021')
# setvar(coreGCP, NA, BETA(), 0.55, "(unitless)")
# setvar(coreGCP, NA, Q10_RH(), 2.2, "(unitless)")
setvar(coreGCP, NA, DIFFUSIVITY(), 1.16, "cm2/s")
setvar(coreGCP,nbp_constrainer$year,"NBP_constrain",nbp_constrainer$nbp,"Pg C/yr")
hector::run(coreGCP,runtodate = 2050)
out <-fetchvars(coreGCP, dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=1)
for (k in c(2, 6)){
core_k <- hector::newcore(ini_file245, name = paste0('245_k',k))
# setvar(coreGCP, NA, BETA(), 0.55, "(unitless)") # have to be consistent with above
# setvar(coreGCP, NA, Q10_RH(),2.2, "(unitless)") # have to be consistent with above
setvar(core_k, NA, DIFFUSIVITY(), 1.16, "cm2/s")
setvar(core_k, nbp_constrainer$year, NBP_CONSTRAIN(), k*nbp_constrainer$nbp, "Pg C/yr")
hector::run(core_k,runtodate = 2050)
fetchvars(core_k,dates=1965:2050,vars=c(GLOBAL_TAS(),
CONCENTRATIONS_CO2(),
NBP(),
LL_OCEAN_UPTAKE(),
HL_OCEAN_UPTAKE())) %>%
mutate(k=k) %>%
bind_rows(out) ->
out
rm(core_k)
}
# Order the variables for faceting
out %>%
mutate(variable = if_else( variable == 'NBP', 'A.NBP',
if_else(variable == 'CO2_concentration', 'B.CO2_concentration',
if_else(variable == 'global_tas', 'C.global_tas',
if_else(variable == 'LL_ocean_uptake', 'D.LL_ocean_uptake',
if_else(variable == 'HL_ocean_uptake', 'E.HL_ocean_uptake',
variable)))))) %>%
mutate(variable = paste0(variable, '~', units)) ->
out2
ggplot(out2 %>% filter(year <= 2000)) +
geom_line(mapping = aes(x=year, y = value, color = k, linetype = scenario)) +
facet_wrap(~variable, scales = 'free_y', nrow=1) +
geom_vline(xintercept = 1994, size = 0.3) +
scale_color_continuous() #+
#ggtitle('diffusivity 1.16 cm2/s, beta=0.55, q10=2.2 \n just changing diff is not enough')
diff: 1.16, beta:.55, q10: 2.2
can’t replicate but ocean uptake in 1993-1995 interesting: in default, both LL and HL ocean are a valley at the same time as NBP valley.
But when that NBP valley is made more extreme, HL and LL become peaks
diff:2.31, beta:0.53, q10:1.76
TODO have plots for both parameterization, not just toggle in code
ALSO TODO: why changing beta and q10 changing plots for same diffusivity (double check side by side) - have to use diff 1.16 to solve with both sets of beta/q10 but identical as should be (think)
is there a rebound when feeding on next time step
ALSO double check gcamlandc - is it passing emissions/land (spikes) or stock/land stock/land might have spikes too, since dividing by land tricky at the changes, esp with lag